REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis 
Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a 
collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

14-06-2006  Technical  Paper 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

Quasi-2D  Unsteady  Flow  Solver  Module  for  Rocket  Engine  and  Propulsion  System 
Simulations  (Preprint) 

5a.  CONTRACT  NUMBER 

FA9300-04-C-0008 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Bryan  T.  Campbell  (Aerojet);  Roger  L.  Davis  (UC  Davis) 

5d.  PROJECT  NUMBER 

502603BQ 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Aerojet 

P.O.  Box  13222 

Sacramento  CA  95813-6000 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL-PR-ED-TP-2  006-189 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/PRS 

5  Pollux  Drive 

Edwards  AFB  CA  93524-7048 

10.  SPONSOR/MONITOR’S 
ACRONYM(S) 

11.  SPONSOR/MONITOR’S 
NUMBER(S) 

AFRL-PR-ED-TP-2  006-189 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  unlimited  (AFRL-ERS-PAS-2006-137) 


13.  SUPPLEMENTARY  NOTES 

Presented  at  the  42nd  AIAA/ASME/SAE/ASEE  Joint  Propulsion  Conference,  Sacramento,  CA,  9-12  July  2006. 


14.  ABSTRACT 

A  new  quasi-two-dimensional  procedure  is  presented  for  the  transient  solution  of  real-fluid  flows  in  lines  and  volumes  including  heat 
transfer  effects.  The  solver  is  targeted  to  the  commercial  dynamic  simulation  software  package  Simulink®  for  integration  into  a  larger  suite  of 
modules  developed  for  simulating  rocket  engines  and  propulsion  systems.  A  Fortran95  code  using  more  conventional  solution  procedures  is 
being  developed  in  parallel  to  provide  verification  test  cases.  The  solution  procedure  for  both  codes  is  coupled  with  a  state-of-the-art  real-fluids 
property  database  so  that  both  compressible  and  incompressible  fluids  may  be  considered  using  the  same  procedure.  The  numerical  techniques 
used  in  this  procedure  are  described.  Test  cases  modeling  transient  flow  of  nitrogen,  water,  and  hydrogen  are  presented  to  demonstrate  the 
capability  of  the  current  technique. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE 

OF  ABSTRACT 

OF  PAGES 

PERSON 

Mr.  Sean  Kenny 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

A 

12 

19b.  TELEPHONE  NUMBER 

(include  area  code) 

N/A 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


Unclassified 


Unclassified 


Unclassified 


Quasi-2D  Unsteady  Flow  Solver  Module  for  Rocket  Engine 
and  Propulsion  System  Simulations  (Preprint) 


Bryan  T.  Campbell 

Aerojet,  Sacramento,  CA,  95813-6000 
and 

Roger  L.  Davis 1 

University  of  California,  Davis,  Davis,  CA,  95616 


A  new  quasi-two-dimensional  procedure  is  presented  for  the  transient  solution  of  real- 
fluid  flows  in  lines  and  volumes  including  heat  transfer  effects.  The  solver  is  targeted  to  the 
commercial  dynamic  simulation  software  package  Simulink®  for  integration  into  a  larger 
suite  of  modules  developed  for  simulating  rocket  engines  and  propulsion  systems.  A 
Fortran95  code  using  more  conventional  solution  procedures  is  being  developed  in  parallel 
to  provide  verification  test  cases.  The  solution  procedure  for  both  codes  is  coupled  with  a 
state-of-the-art  real-fluids  property  database  so  that  both  compressible  and  incompressible 
fluids  may  be  considered  using  the  same  procedure.  The  numerical  techniques  used  in  this 
procedure  are  described.  Test  cases  modeling  transient  flow  of  nitrogen,  water,  and 
hydrogen  are  presented  to  demonstrate  the  capability  of  the  current  technique. 
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Nomenclature 

=  cross-sectional  area 
=  perimeter 

=  relative  acceleration  between  inertial  and  absolute  frames  of  reference 
=  fluid  speed  of  sound 
=  wall  material  specific  heat 
=  flow  passage  hydraulic  diameter 
=  fluid  total  energy  =  p  (e  +  if/2  +  gz) 

=  fluid  internal  energy 
=  fluid  friction  coefficient 

=  vector  of  fluid  mass,  momentum,  and  energy  fluxes 
=  gravitational  acceleration 
=  forced  convection  heat  transfer  coefficient 

=  fluid  static  enthalpy  referenced  to  fluid  normal  boiling  point  at  latm 
=  cell  index 
=  minor  loss  coefficient 
=  fluid  thermal  conductivity 
=  wall  material  thermal  conductivity 
=  domain  length 

=  Nusselt  number  based  on  hydraulic  diameter  (Nud  -  hD/k) 

=  fluid  static  pressure 

=  fluid  Prandtl  number  evaluated  at  film  temperature  (Pr  =  cp/i  / k) 

=  heat  flux  from  wall  to  fluid 

=  local  flow  Reynolds  number  based  on  hydraulic  diameter  (ReD  =  pnD//j) 

=  local  flow  Reynolds  number  based  on  hydraulic  diameter  evaluated  at  film  temperature 
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=  source  vector 
=  time 

=  fluid  bulk  temperature 
=  film  temperature  (7}  =  (Tw+TB)/2) 

=  wall  temperature 

=  fluid  state  vector 

=  fluid  absolute  velocity 

=  relative  velocity  of  moving  frame  of  reference 

=  fluid  cell  volume 

=  distance  along  solution  domain 

=  elevation  above  reference  plane 

=  length  of  computational  cell  i 

=  wall  roughness  height 

=  fluid  density 

=  wall  material  density 

=  shear  force  between  wall  and  fluid  (friction) 

=  fluid  molecular  viscosity 


I.  Introduction 

THIS  paper  addresses  development  of  a  quasi-two  dimensional  (2D),  unsteady,  two-phase  flow  solver  with  real- 
fluid  properties  suitable  for  use  in  system- level  transient  simulations  of  hydraulic  and  pneumatic  systems  and,  in 
particular,  rocket  engine  and  propulsion  systems.  The  flow  solver  module  described  here  is  an  integral  part  of  a 
rocket  engine  system  modeling  tool  currently  under  development.  The  solver  is  geared  towards  modeling  the 
dynamic  behavior  of  fluid-filled  lines  and  passages  (i.e.,  the  solution  domain  is  much  larger  in  one  dimension  than  in 
the  others)  where  heat  transfer  and  phase  change  are  important.  Existing  software  models  available  for  such 
simulations  are  typically  based  on  lumped  analysis  approaches,  which  do  not  fully  capture  the  flow  dynamics, 
mainly  due  to  the  elimination  of  the  unsteady  momentum  terms.  Although  so-called  “continuity”  waves  can  be 
captured  using  a  lumped  analysis,  neglecting  these  terms  leads  to  an  inability  to  capture  true  “dynamic”  waves 
required  for  simulating  such  phenomenon  as  water-hammer  and  pressure  surge1.  Since  these  phenomena  play  an 
important  role  in  the  operation  and  testing  of  liquid  rocket  engines  and  propulsion  systems,  a  method  that  captures 
these  effects  is  desired. 

Several  approaches  for  simulating  the  dynamic  behavior  of  such  “fluid-transmission”  lines  have  been  reported. 
The  “lumped-analysis”  approach  treats  a  flow  passage  as  a  series  of  fluid  control  volumes  that  conserve  mass  and 
energy  linked  by  flow  resistance  elements  that  compute  the  flow  between  the  volumes2.  While  this  approach  does 
conserve  momentum  in  a  quasi-steady  sense  at  the  flow  resistances,  the  unsteady  momentum  term  in  the  governing 
equations  is  neglected.  Another  approach  uses  the  Method  of  Characteristics3,  which  is  a  general  method  for  solving 
hyperbolic  partial  differential  equations.  The  governing  equations  for  fluid  flow  are  compatible  with  this  method  and 
it  has  been  used  for  simulation  of  fluid  transmission  lines4.  While  the  unsteady  momentum  terms  are  retained  using 
this  method,  other  problems,  particularly  at  the  boundaries  of  components,  make  it  problematic  to  apply  to  a 
modular  system-level  simulation  tool.  Modal  methods  have  also  been  used  when  solution  in  the  frequency  domain  is 
possible5.  This  technique  represents  the  pressure  and  flow  distribution  in  the  flow  domain  as  a  sum  of  an  infinite 
series  of  mode  shapes,  similar  to  a  Fourier  series  solution.  While  this  method  presents  an  elegant  and  efficient 
method  for  simulating  simple  flow  situations  (e.g.,  incompressible  flow,  inviscid  flow,  laminar  flow,  etc.),  the 
addition  of  turbulent  flow,  real-fluid  properties,  heat  transfer,  and  phase  change  complicate  application  of  the 
method  and  reduce  its  attractiveness. 

The  current  technique  focuses  on  application  of  standard  finite -volume  methods  to  the  development  of  an 
unsteady  flow  solver  capable  of  capturing  wave  dynamics,  real-fluid  effects,  heat  transfer,  and  phase  change.  In 
addition,  the  solution  method  must  be  computationally  efficient  enough  to  operate  as  a  module  in  a  larger  system- 
level  simulation.  The  following  sections  describe  the  modeling  approach,  numerical  methodologies,  and  test  cases 
that  have  been  used  during  development  of  this  model. 


II.  Approach 

The  intention  of  the  model  developed  here  is  to  represent  fluid  lines  and  flow  passages  as  a  quasi-2D  domain. 
The  intended  components  to  be  modeled  using  this  method  are  fluid  lines  and  internal  flow  passages  where  the 
length  of  the  domain  is  much  larger  than  the  hydraulic  diameter  of  the  domain.  Flow  separations  and  non-axial 
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velocities  are  minimal  for  these  types  of  components;  hence,  the  quasi-2D  assumption  is  valid.  The  solver  is  targeted 
to  the  commercial  Simulink8  dynamic  simulation  software  package  from  The  MathWorks  for  integration  into  a 
larger  suite  of  modules  developed  for  simulating  rocket  engines  and  propulsion  systems.  Simulink8  was  selected 
since  it  offers  a  wide  range  of  capabilities,  over  a  dozen  robust  differential  equation  solvers,  extensive 
documentation  and  technical  support,  a  modern  graphically  based  modeling  paradigm,  an  existing  user  community 
across  many  disciplines,  and  commercially  funded  code  development  and  maintenance.  A  Fortran95  code  is  also 
being  developed  in  parallel'1  to  provide  test  cases  for  the  Simulink8  module. 

The  solver  is  being  developed  to  model  the  effects  of  variable  area  along  the  length  of  the  domain,  friction, 
minor  losses  (e.g.,  bends),  real-fluid  and  two-phase  flow  effects,  gravity  and  acceleration  effects,  heat  transfer,  and 
the  capability  to  produce  unsteady  and  steady-state  solutions.  Fluid  equation  of  state  closure  is  provided  by  the 
REFPROP  fluid  property  database7  available  from  the  National  Institute  for  Standards  and  Testing  (NIST).  This 
database  uses  state-of-the-art  Equation  of  State  models  to  fully  describe  “real-fluid”  properties  over  a  wide  range  of 
thermodynamic  conditions,  including  liquid,  vapor,  mixed  phase,  and  supercritical  fluid  regimes.  Properties  that 
completely  define  the  fluid  thermodynamic  state,  as  well  as  transport  properties,  are  available  as  a  function  of  any 
two  thermodynamic  parameters.  Validated  fluid  models  for  over  80  pure  fluids  and  over  180  fluid  mixtures  are 
available  in  the  database.  The  database  is  accessed  through  a  suite  of  Fortran  subroutines  that  were  compiled  into  the 
codes  and  obtain  fluid  equation  of  state  model  parameters  from  external  files. 

Friction  and  heat  transfer  are  modeled  as  source  terms  in  the  fluid  governing  equations.  This  approach  allows  the 
flow  to  be  modeled  as  one-dimensional  (ID)  and  facilitates  computational  efficiency.  Friction  (i.e.,  viscous  losses 
and  losses  due  to  bends  and  fittings)  and  heat  transfer  coefficients  are  obtained  from  suitable  correlations  between 
the  flow  variables  and  source  terms. 


III.  Governing  Equations 

The  modeling  procedure  solves  differential  equations  representing  the  physics  of  the  problem.  Equations  are 
presented  below  for  modeling  the  fluid  dynamics,  wall  temperature,  and  source  terms  coupling  the  two. 


A.  Fluid 

The  differential  equations  governing  the  dynamics  of  the  fluid  consist  of  the  quasi  two-dimensional  continuity, 
Navier-Stokes  and  energy  equations  simplified  to  model  viscous  effects  and  general  acceleration  in  a  non-inertial 
frame  using  source  terms8. 
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B.  Wall 

The  wall  bounding  the  fluid  computational  domain  is  modeled  using  standard  control  volume  methods  for  heat 
transfer9.  Forced  convection  heat  transfer  between  the  fluid  and  wall  as  well  as  heat  conduction  within  the  wall  are 
included  in  the  model.  Fleat  conduction  between  the  wall  and  fluid  is  neglected  due  to  its  small  magnitude  compared 
to  forced  convection  heat  transfer.  In  addition,  conduction  in  the  wall  is  modeled  as  1-D  neglecting  radial 
conduction.  These  effects  can  be  added  to  the  model  as  required. 

Based  on  these  assumptions,  the  governing  equation  for  heat  transfer  within  the  wall  is9: 


dt 
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Eq.  (3)  does  not  include  the  contribution  due  to  forced  convection  heat  transfer  between  the  fluid  and  the  wall.  It  is 
more  convenient  to  add  this  term  below  to  the  discretized  equations. 

C.  Source  Terms 

Equations  (1)  through  (3)  above  require  several  source  terms  that  are  a  function  of  the  wall  and  fluid  conditions. 
The  wall  shear  stress  can  be  related  to  a  suitable  friction  factor,  f  which  is  a  function  of  the  local  Reynolds  number 
and  wall  surface  roughness.  For  single-phase  flow  in  pipes,  the  Churchill  correlation10  can  be  used  to  determine  wall 
friction  factors  using  the  equation: 
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Similarly,  forced  convection  heat  transfer  between  the  wall  and  fluid  can  be  obtained  using  a  suitable  correlation. 
For  single-phase  flow  in  pipes,  the  Colburn  correlation"  can  be  used  to  determine  the  forced  convection  heat 
transfer  from  the  wall  to  fluid: 


qw  =  h  ■  Ap  ■  {Tw  - Tb),Nud=  —  =  0.023  •  Re0/  ■  Prf  (5) 

The  source  term  correlations  given  here  are  suitable  for  single-phase  flows  of  non-specific  fluids.  Other  fluid- 
specific  and/or  two-phase  correlations  can  be  substituted  to  model  other  types  of  flows  as  required. 

D.  Boundary  Conditions 

Boundary  conditions  are  applied  each  time-step  at  the  inlet  and  exit  of  each  line  in  a  network.  These  boundary 
conditions  are  specified  depending  on  network  connectivity  and  consist  of  a  global  inlet,  global  exit,  or  junction 
boundary  condition.  At  the  global  inlet(s)  to  the  network,  the  stagnation  enthalpy  and  entropy  are  specified  based  on 
a  given  stagnation  pressure  and  stagnation  temperature  and  the  predicted  velocity.  Resulting  fluid  density  and 
internal  energy  are  determined  from  REFPROP,  thus  allowing  for  determination  of  the  fluid  states  at  the  inlet.  At  the 
global  exit(s)  to  the  network,  the  static  pressure  is  specified.  The  internal  energy  and  velocity  states  along  with  this 
specified  pressure  are  used  with  REFPROP  to  determine  density  and  total  energy.  The  prescribed  inlet  and  exit 
conditions  may  be  specified  as  functions  of  time.  At  line  junctions,  the  contributions  of  the  time-rate  changes  in  the 
primary  variables  contained  in  U  of  Eq.  (1)  and  Tw  of  Eq.  (3)  are  summed  to  give  total  time -rate  change  at  the 
junction  node.  For  each  junction  of  a  given  line,  the  neighboring  line  numbers  and  boundaries  (i.e.,  inlet  or  exit)  are 
stored  to  make  implementation  of  this  boundary  condition  straightforward. 

IV.  Numerical  Methods 

Equations  (1)  through  (3)  may  be  solved  with  a  variety  of  numerical  schemes  that  are  second-order  accurate 
in  space  and  time  for  “steady”  or  transient  flow.  The  flow  solver  is  implemented  using  an  ANSI  C  S-Function, 
which  is  a  specialized  program  for  use  in  Simulink6  models  that  contains  specific  subroutines  required  during  model 
execution.  The  model  states  are  the  flow  and  wall  variables  contained  in  the  U  vector  and  Tw  at  each  computational 
node  and  the  main  subroutine  of  the  S-Function  computes  the  time  derivatives  of  the  flow  states  ( dU/dt  and  dTJdt). 
Once  time  derivatives  of  the  system  states  are  computed,  standard  Simulink1'  solvers  are  used  to  integrate  the 
derivatives  and  compute  the  time  history  of  the  model  states.  Multiple  numerical  schemes  have  been  implemented  in 
the  Fortran95  code  to  provide  flexibility  and  robustness  when  solving  for  different  flow  conditions.  This  version  of 
the  code  uses  solution  methods  similar  to  those  discussed  here  and  is  described  in  more  detail  elsewhere6. 

A.  Computational  Grid 

The  governing  differential  equations  are  solved  on  a  1-D  grid  with  varying  area  as  shown  in  Fig.  1.  In  general, 
the  computational  domain  has  N  computational  cells.  The  fluid  states  ( U)  are  located  at  the  cell  faces,  as  shown  in 
the  figure.  The  wall  states  (Tw)  are  located  in  the  wall  at  the  midpoint  of  the  cells  as  shown.  This  staggered 
arrangement  was  selected  to  simplify  application  of  boundary  conditions  for  the  fluid  and  wall.  Thus,  for  N  cells, 
there  are  3(N+1)  fluid  states  and  N  wall  states. 
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Wall  Domain 
Fluid  Domain 


Figure  1.  Computational  grid. 


B.  Discretization 

As  discussed  above,  the  main  subroutine  of  the  flow  solver  S-Function  computes  the  model  state  derivatives 
( dU/dt  and  dTJdi)  at  each  computational  node.  A  finite -volume  method  is  used  to  compute  the  fluid-state  fluxes  at 
the  center  of  each  fluid  cell,  so  for  cell  1+V2  in  Fig.  1: 
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The  fluid-state  fluxes  from  each  cell  are  distributed  to  the  fluid-state  derivatives  computed  at  the  neighboring  faces: 
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A  similar  method  is  used  to  compute  wall-state  fluxes  at  the  center  of  each  wall  cell9;  however,  since  the  wall 
states  are  at  the  center  of  the  wall  cells — no  distribution  of  fluxes  is  required: 
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where  the  wall  material  specific  heat  (c„,)  is  a  function  of  local  wall  temperature,  and  k+  and  k.  are  weighted  wall 
thermal  conductivity  values  for  neighboring  cells  computed  using  the  harmonic  mean  discussed  in  Reference  9. 
Note  that  the  forced  convection  heat  transfer  between  the  wall  and  fluid  has  been  introduced. 

C.  Temporal  Solution 

Once  the  state  derivatives  are  computed  at  each  computational  node,  standard  SimulinkR  solvers  are  used  to 
solve  the  discretized  equations  as  a  function  of  time.  Four  variable-step  solvers  were  investigated  for  solving 
Eqs.  (1)  through  (3)  on  the  computational  domain,  all  of  which  are  suited  to  stiff  systems  of  differential  equations 
like  those  considered  here.  The  most  rigorous  of  these  is  the  implicit  solver  ode23s,  which  is  a  one-step  solver  based 
on  a  modified  Rosenbrock  formula  of  order  two12.  Odel5s  is  a  multi-step  solver  based  on  Numerical 
Differentiation  Formulas  (NDFs)12,  which  are  related  to  the  backward  differentiation  formulas  (also  known  as  Gear's 
method)  but  are  more  efficient13.  Ode23t  is  an  implementation  of  the  trapezoidal  rule  using  a  "free"  interpolant  that 
is  effective  for  moderately  stiff  problems13.  Ode23tb  uses  the  TR-BDF214  two-stage  procedure  where  the  first  stage 
is  a  trapezoidal  rule  step  and  the  second  stage  is  a  backward  differentiation  formula  of  order  two13.  References  12 
and  13  provide  additional  details. 
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V.  Results 


Several  test  cases  have  been  generated  to  evaluate  the  merits  of  the  procedure  described  here.  Solutions  were 
obtained  with  both  the  Simulink®  code  and  Fortran95  version  of  the  solver.  Results  are  presented  below. 

A.  Compressible  Flow  in  a  Pipe 

A  transient  pipe  flow  simulation  of  gaseous  nitrogen  in  a  pipe  has  been  performed  to  demonstrate  the  capability 
of  the  present  procedure  for  compressible  flows.  A  0.254  m  (10  in.)  pipe  with  a  0.0254  m  (1  in.)  diameter  and  90- 
deg  elbow  at  0.127  m  (5  in.)  was  used  for  these  simulations  as  shown  in  Fig.  2.  The  surface  roughness  of  the  pipe 
was  2.54xl0"6  m  (0.0001  in.),  a  typical  value  for  commercially  smooth  pipe.  The  inlet  total  as  well  as  initial  exit 
static  pressure  were  held  at  345  kPa  (50  psia).  The  inlet  total  temperature  was  held  at  300°K  (540°R).  Thus,  the 
initial  velocity  of  flow  in  the  pipe  was  zero.  At  0.001  seconds,  the  exit  static  pressure  was  reduced  to  331  kPa  (48 
psia)  over  0.001  seconds  (linear  ramp).  This  reduction  in  pressure  simulates  a  fast-acting  valve  opening  downstream 
of  the  line,  which  initiates  a  series  of  expansion  and  compression  waves  in  the  pipe  and  an  increase  in  flow  velocity. 
At  steady-state,  static  pressure  in  the  pipe  settles  to  a  non-uniform  distribution  along  the  length  of  the  pipe  with  inlet 
pressure  higher  than  exit  due  to  friction  losses.  The  flow  velocity  grows  asymptotically  until  it  reaches  a  steady-state 
distribution.  The  simulation  was  run  to  steady-state  when  both  pressure  and  velocity  distributions  remain  lime- 
independent.  A  total  of  32  cells  with  uniform  spacing  were  used  along  the  pipe  to  discretize  the  governing  equations. 


Figure  2.  Simple  pipe  geometry. 

Figure  3  shows  the  predicted  pressure,  velocity,  and  enthalpy  as  a  function  of  time  using  the  implicit  Simulink'' 
solver  ode23s.  Sinusoidal  variation  of  pressure  with  time  is  evidence  of  the  expansion  and  compression  waves 
propagating  through  the  pipe.  Time-variation  of  the  velocity  corresponds  to  these  waves.  Velocity  increases  until  it 
reaches  an  average  value  of  65.9  m/s  (216.3  ft/s).  Figure  4  shows  the  pressure,  velocity,  and  enthalpy  distributions  at 
steady-state.  Pressure  and  enthalpy  decrease  along  the  length  of  the  pipe  due  to  viscous  effects,  with  a  large  jump  at 
0.127  m  (5  in.)  due  to  the  minor  loss  in  the  elbow.  Figure  5(a)  compares  transient  pressure  predictions  of  the  model 
at  the  mid-point  of  the  pipe  using  four  different  Simulink1'  solvers.  The  solvers,  based  on  the  NDF  equations 
(odel5s),  trapezoidal  rule  (ode23t),  and  TR-BDF2  solver  (ode23tb),  give  very  similar  results  to  the  implicit  solver 
(ode23s).  However,  these  solvers  run  in  about  4  percent  of  the  computational  time  required  for  the  implicit  solver 
(7,443  seconds  for  the  0.02  second  simulation),  making  them  attractive  for  quick  turn-around  analyses.  The  NDF 
solver  only  slightly  under-predicts  the  pressure  oscillation  frequency  after  three  cycles  while  the  other  solvers  show 
better  agreement  to  frequency  and  amplitude  of  the  implicit  solver  results.  Finally,  Fig.  5(b)  compares  the  pressure 
solution  obtained  with  the  implicit  solver  (ode23s)  with  the  solution  obtained  using  the  explicit  version  of  the 
Fortran  code  discussed  in  Reference  6.  The  two  solvers  produce  very  similar  results  with  only  small  differences 
between  pressure  wave  amplitude  and  frequency.  The  0.02  second  Fortran  simulation  runs  in  about  135  seconds 
using  the  explicit  solver,  or  in  about  half  the  time  of  the  faster  Simulink  solvers.  This  comparison  gives  an  estimate 
of  the  computational  overhead  associated  with  using  Simulink  for  this  type  of  calculation. 
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Figure  3.  Nitrogen  transient  pipe  flow. 
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Figure  4.  Nitrogen  pipe  flow  at  steady-state. 


(a)  Four  Simulink  Solvers 


(b)  Comparison  of  ODE23s  to  Explicit  FORTRAN  Solver 
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Figure  5.  Transient  solution  with  various  solvers. 

A  hand  calculation1’  for  the  steady-state  case  assuming  incompressible  flow  yields  a  static  pressure  drop  of 
5.4095  kPa  (0.7846  psid)  and  flowrate  of  0.1275  kg/s  (0.2810  lbm/s).  The  simulation  yields  a  steady-state  static 
pressure  drop  of  5.5480  kPa  (0.8047  psid)  and  flowrate  of  0.1264  kg/s  (0.2788  lbm/s).  Thus,  the  simulation  results 
differ  from  the  hand  calculations  by  2.6  percent  in  pressure  and  0.8  percent  in  flow.  However,  the  hand  calculation  is 
expected  to  have  some  error  due  to  the  incompressibility  assumption.  Thus,  the  steady-state  simulation  results  are 
verified  by  independent  hand  calculations. 

The  period  of  pressure  oscillations  can  be  compared  to  the  period  expected  for  frictionless  flow  (e.g.,  organ  pipe 
modes)  as  a  check  of  the  transient  response  of  the  computational  solution.  The  pressure  waves  inside  the  pipe  travel 
at  slightly  less  than  the  speed  of  sound,  where  the  difference  is  a  result  of  friction  losses.  Thus,  pressure  oscillations 
should  occur  close  to  the  natural  frequency  of  the  pipe  (c/2L),  or  696Hz  for  this  case.  The  oscillation  period  of  the 
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first  wave  in  the  simulation  is  0.00147  seconds,  which  is  slightly  higher  than  the  period  corresponding  to  the  natural 
frequency  of  the  pipe  (0.00144  seconds).  The  period  of  subsequent  waves  is  longer  due  to  the  increasing  effects  of 
friction  losses  as  flow  velocity  increases.  Thus,  the  method  produces  transient  results  in  keeping  with  physical 
expectations. 


B.  Incompressible  Flow  in  a  Pipe 

The  same  configuration  described  above  in  Fig.  2  was  used  for  water  flow  simulation  to  demonstrate  the 
capability  to  predict  transient  flows  in  essentially  incompressible  fluids.  The  critical  difference  for  incompressible 
fluids  is  that  wave  speeds  are  much  larger  than  in  compressible  fluids,  resulting  in  smaller  time  steps  for  numerical 
stability.  The  same  time -dependent  boundary  conditions  as  discussed  above  are  imposed  on  the  inlet  and  outlet 
planes  of  the  pipe.  Figures  6  and  7  show  the  transient  and  steady-state  solutions  for  the  water  flow  simulation 
computed  using  solver  ode23tb  (the  TR-BDF2  solver)  with  a  maximum  time  step  of  1  x  10  s  seconds.  Small 
amplitude  pressure  oscillations  are  evident  throughout  the  duration  of  the  simulation,  and  the  flow  relaxes  gradually 
to  its  equilibrium  state.  The  time  to  reach  steady-state  for  water  is  much  longer  than  that  required  for  nitrogen  due  to 
the  incompressibility  of  the  fluid.  At  steady-state,  the  flow  velocity  is  nearly  constant,  as  expected  for  a  nearly 
incompressible  fluid.  Figure  8  shows  the  first  0.004  seconds  of  static  pressure  and  enthalpy  variation  (velocity 
variation  is  initially  negligible)  and  demonstrates  good  agreement  of  the  Simulink  results  to  the  explicit  version  of 
the  Fortran  code  for  oscillation  frequency  and  magnitude. 

As  with  the  nitrogen  case,  the  steady-state  water  solution  can  be  compared  to  a  hand  calculation  for  the  same 
conditions.  Both  the  hand  calculation  and  simulation  yield  flowrates  of  2.0508  kg/s  (4.5213  lbm/s)  and  static 
pressure  drops  of  5.5762  kPa  (0.8088  psid).  The  agreement  between  the  hand  calculation  and  simulation  is  better  for 
water  since  the  assumption  of  incompressible  flow  is  more  correct  for  water  compared  to  gaseous  nitrogen. 


(a)  Static  Pressure  (b)  Velocity  (c)  Static  Enthalpy 
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Figure  6.  Water  transient  pipe  flow. 
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Figure  7.  Water  pipe  flow  at  steady-state. 
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Solid  =  Simulink 
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Figure  8.  Comparison  of  water  transient  results  with  Fortran  code. 

C.  Cryogenic  Fluids  With  Heat  Transfer 

The  current  procedure  can  also  be  used  to  solve  for  transient  behavior  of  cryogenic  fluids,  including  heat  transfer 
between  fluid  and  pipe.  A  simple  case  similar  to  the  one  described  above  was  modeled  with  liquid  parahydrogen  to 
demonstrate  this  capability.  A  0.6096  m  (24  in.)  long  straight  pipe  with  an  inner  diameter  of  0.0254  m  (1  in.)  and 
surface  roughness  of  2.54xl0'6  m  (0.0001  in.)  was  used.  The  initial  fluid  temperature  was  set  at  22.22  °K  (40°  R) 
and  initial  pressure  was  345  kPa  (50  psia).  Inlet  total  pressure  was  held  at  345  kPa  (50  psia)  and  inlet  total 
temperature  was  held  at  22.22  °K  (40  °R).  Thus,  the  initial  flow  velocity  in  the  pipe  was  zero.  The  copper  pipe  had  a 
0.003175  m  (0.125  in.)  thick  wall,  which  was  initially  set  to  the  same  temperature  as  the  fluid.  At  0.001  seconds,  the 
exit  static  pressure  was  reduced  to  331  kPa  (48  psia)  over  0.001  seconds  (linear  ramp).  The  pipe  wall  was  subjected 
to  zero  heat  transfer  at  the  inlet  and  outlet  boundaries.  The  conditions  for  this  problem  were  chosen  to  avoid  two- 
phase  flow  since  this  capability  is  currently  not  implemented.  The  problem  was  modeled  using  the  TR-BDF2  solver 
(ode23tb)  with  a  maximum  time  step  of  lxlO"5  seconds  and  32  computational  cells  along  the  domain.  While 
temperature  change  of  the  wall  is  expected  to  be  small  for  this  case,  it  demonstrates  the  capability  to  model 
conjugate  heat  transfer  using  this  procedure. 

Figures  9  and  10  show  the  transient  and  steady-state  fluid  solutions  for  the  liquid  hydrogen  flow  simulation. 
Pressure  oscillations  of  moderate  amplitude  are  evident  throughout  the  duration  of  the  simulation,  and  the  flow 
relaxes  gradually  to  its  equilibrium  state.  Since  liquid  hydrogen  compressibility  is  lower  than  gaseous  nitrogen  and 
higher  than  liquid  water,  its  behavior  is  between  the  previous  cases.  The  time  to  reach  steady-state  is  similar  to  the 
water  case.  At  steady-state,  the  flow  velocity  is  nearly  constant,  as  expected  for  a  fluid  with  low  compressibility. 
Figure  11  shows  the  first  0.006  seconds  of  static  pressure  and  enthalpy  variation  (velocity  variation  is  initially 
negligible)  and  demonstrates  good  agreement  of  the  Simulink  results  to  the  explicit  version  of  the  Fortran  code  for 
oscillation  frequency  and  magnitude. 


(a)  Static  Pressure  (b)  Velocity  (c)  Static  Enthalpy 
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Figure  10.  Liquid  hydrogen  pipe  flow  at  steady-state. 


Figure  11.  Comparison  of  hydrogen  transient  results  with  Fortran  code. 

Figures  12  and  13  show  the  pipe  wall  temperature  history  and  distribution  during  the  simulation.  As  shown  in 
Fig.  12,  wall  temperature  decays  as  flow  accelerates  and  static  temperature  of  the  liquid  hydrogen  is  reduced.  Wall 
temperature  begins  to  change  first  at  the  pipe  outlet  since  the  fluid  velocity  increases  there  first.  The  wall 
temperature  reaches  a  steady-state  distribution  within  0.3  seconds.  Figure  13  shows  the  wall  temperature  distribution 
along  the  length  of  the  pipe  at  several  times  during  the  simulation.  Initially  (t — 0),  the  pipe  wall  is  at  a  constant 
temperature.  As  flow  is  established  at  the  exit  of  the  pipe,  wall  temperature  drops  due  to  forced  convection  heat 
transfer.  As  flow  velocity  increases  during  the  simulation,  wall  temperature  continues  to  drop.  At  end  of  the 
simulation  (t=0.3  seconds),  the  pipe  wall  has  reached  a  temperature  distribution  in  equilibrium  with  the  fluid  static 
temperature  distribution,  which  is  shown  as  circular  symbols  in  the  plot. 


Figure  13.  Wall  temperature  distribution. 
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As  with  the  previous  cases,  the  steady-state  solution  can  be  compared  to  a  hand  calculation  for  the  same 
conditions.  The  hand  calculation  yields  a  flowrate  of  0.6109  kg/s  (1.3467  lbm/s)  and  static  pressure  drop  of 
3.2222  kPa  (0.4673  psid).  At  steady-state,  the  simulation  indicates  a  flowrate  of  0.6108  kg/s  (1.3467  lbm/s)  and 
static  pressure  drop  of  3.2230  kPa  (0.4675  psid),  or  percent  differences  of  0.004  and  0.02,  respectively.  The 
agreement  between  the  incompressible  hand  calculation  and  simulation  is  better  for  liquid  hydrogen  compared  to 
gaseous  nitrogen,  but  not  as  good  as  water  since  the  compressibility  of  liquid  hydrogen  falls  between  liquid  water 
and  gaseous  nitrogen. 


VI.  Conclusion 

A  new  baseline  procedure  has  been  developed  for  the  numerical  solution  of  transient  quasi-2D  flow  in  lines  and 
volumes  including  the  effects  of  friction,  minor  losses,  real-fluid  properties,  and  heat  transfer  between  the  fluid  and 
bounding  wall.  The  solver  is  specifically  targeted  to  rocket  engine  and  propulsion  system  modeling,  but  can  be 
applied  to  other  types  of  systems.  The  procedure  has  been  implemented  in  both  Matlab/Simulink1'  and  Fortran95 
using  a  range  of  explicit  and  implicit  numerical  schemes  for  the  solution  of  the  differential  equations  to  provide 
flexibility.  Real-fluid  properties  are  obtained  from  a  NIST  fluid  property  database,  allowing  both  compressible  and 
incompressible  flow  problems  to  be  treated  using  the  same  procedure.  For  cases  with  wall  heat  transfer,  the  heat 
conduction  equation  is  solved  simultaneously  with  the  fluid  equations  to  provide  conjugate  solutions.  Various 
compressible  and  incompressible  unsteady  flow  test  cases  are  presented  to  verify  solution  accuracy  against 
analytical  solutions.  Future  work  will  focus  on  extension  of  the  procedure  to  handle  variable  flow  area,  two-phase 
flow,  and  acceleration;  increasing  computational  efficiency  through  the  use  of  parallel  computing;  and  application  of 
the  module  for  system-level  analyses. 
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